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ABSTRACT 

Photometry of stars from the K2 extension of NASA’s Kepler mission is 
afflicted by systematic effects cansed by small (few-pixel) drifts in the telescope 
pointing and other spacecraft issues. We present a method for searching K2 light 
curves for evidence of exoplanets by simultaneously htting for these systematics 
and the transit signals of interest. This method is more computationally expen¬ 
sive than standard search algorithms but we demonstrate that it can be efficiently 
implemented and used to discover transit signals. We apply this method to the 
full Campaign 1 dataset and report a list of 36 planet candidates transiting 31 
stars, along with an analysis of the pipeline performance and detection efficiency 
based on artihcial signal injections and recoveries. For all planet candidates, we 
present posterior distributions on the properties of each system based strictly on 
the transit observables. 
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1. Introduction 

The Kepler Mission was incredibly snccessfnl at finding transiting exoplanets in the 
light curves of stars. The Mission has demonstrated that it is possible to routinely measure 
signals in stellar light curves at the part-in-10^ level. Results from the primary mission 
include the detection of planet transits with depths as small as 12 parts per million (Barclay 
et al. 2013). 

The noise floor for Kepler data is often quoted as 15 parts per million (ppm) per six 
hours of observations (Gilliland et al. 2011). Although they generally do not interfere with 
searches for transiting planets, larger systematic effects exist on different timescales. One of 
the most serious of these is spacecraft pointing: If the detector fiat-field is not known with 
very high accuracy, then tiny changes to the relative illumination of pixels caused by a star’s 
motion in the focal plane will lead to changes in the measured or inferred brightness of the 
star. 


The great stability of the original Kepler Mission came to an end with the failure of a 
critical reaction wheel. The K2 Mission (Howell et al. 2014) is a follow-on to the primary 
Mission, observing about a dozen fields near the ecliptic plane, each for ~ 75 days at a time. 
Because of the degraded spacecraft orientation systems, the new K2 data exhibit far greater 
pointing variations—and substantially more pointing-induced variations in photometry— 
than the original Kepler Mission data. This makes good data-analysis techniques even more 
valuable. 

Good photometry relies on either a near-perfect fiat-field and pointing model or else 
data-analysis techniques that are insensitive to these instrument properties. The fiat-field 
for Kepler was measured on the ground before the launch of the spacecraft, but is not 
nearly as accurate as required to make pointing-insensitive photometric measurements at 
the relevant level of precision. In principle direct inference of the fiat-field might be possible; 
however, because point sources are observed with relatively limited spacecraft motion, and 
only a few percent of the data are actually stored and downloaded to Earth, there isn’t 
enough information in the data to derive or infer a complete or accurate fiat-field map. 
Therefore, work on K2 is sensibly focused on building data-analysis techniques that are 
pointing-insensitive. 

Previous projects have developed methods to work with K2 data. Both Vanderburg & 
Johnson (2014) and Armstrong et al. (2014) extract aperture photometry from the pixel data 
and decorrelate with image centroid position, producing light curves for each star that are 
“corrected” for the spacecraft motion. These data have produced the first confirmed planet 
found with K2 (Vanderburg et al. 2014). Both Aigrain et al. (2015) and Grossfield et al. 
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(2015) use a Gaussian Process model for the measured flux, with pointing measurements as 
the inputs, and then “de-trend” using the mean prediction from that model. Other data- 
driven approaches have been developed and applied to the data from space missions (for 
example, Ofir et ah 2010; Stumpe et ah 2012; Smith et ah 2012; Petigura et ah 2013a; Wang 
et ah 2015) and ground-based surveys (for example, Kovacs et ah 2005; Tamuz et ah 2005; 
Berta et ah 2012) but they have yet to be generalized to K2. 

In all of these light-curve processing methodologies, the authors follow a traditional 
procedure of “correcting” or “de-trending” the light curve to remove systematic and stellar 
variability as a step that happens before the search for transiting planets. Fit-and-subtract 
is dangerous: Small signals, such as planet transits, can be partially absorbed into the best- 
£t stellar variability or systematics models, making each individual transit event appear 
shallower. In other words, the traditional methods are prone to over-fitting. Because over¬ 
fitting will in general reduce the amplitude of true exoplanet signals, small planets that 
ought to appear just above any specific signal-to-noise or depth threshold could be missed 
because of the de-trending. This becomes especially important as the amplitude of the noise 
increases. 

The alternative to this approach is to simultaneously fit both the systematics and the 
transit signals. Simultaneous htting can push the detection limits to lower signal-to-noise 
while robustly accounting for uncertainties about the systematic trends. In particular, it 
permits us to marginalize over choices in the noise model and propagate any uncertainties 
about the systematic effects to our confidence in the detection. This marginalization ensures 
that any conclusions we come to about the exoplanet properties are conservative, given the 
freedom of the systematics model. 

In this Article we present a data-analysis technique for exoplanet search and character¬ 
ization that is insensitive to spacecraft-induced trends in the light curves. We assume that 
the dominant trends in the observed light curves in each star are caused by the spacecraft 
and are, therefore, shared with other stars. We reduce the dimensionality by running PGA 
on stellar light curves to obtain the dominant modes. The search for planets proceeds by 
modeling the data as a linear combination of 150 of these basis vectors and a transit model. 
Our method builds on the ideas behind previous data-driven de-trending procedures such 
as the Kepler pipeline pre-search data conditioning (PDC; Stumpe et ah 2012; Smith et ah 
2012), but (because of our simultaneous fitting approach) we can use a much more flexible 
systematics model while being less prone to over-fitting. 

The methods developed within this paper are highly relevant to both K2 and the up¬ 
coming TESS mission (Ricker et ah 2014). TESS will feature pointing precision of ~ 3 arc- 
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seconds^, similar to the level of pointing drift with K2. Moreover, the typical star will be 
only observed for one month at a time, and the typical transit detection will be at a similar 
signal-to-noise ratio as with K2. 

Catalogs of transiting planets found in the K2 data will be important to better under¬ 
stand the physical properties, formation, and evolution of planetary systems. These planets, 
especially when they orbit bright or late-type stars, will be useful targets for ground-based 
and space-based follow-up, both for current facilities and those planned in the near future 
such as JWST. They will also deliver input data for next-generation population inferences 
(Foreman-Mackey et ah 2014), especially for the population of planets around cool stars (for 
example. Dressing & Charbonneau 2015). 

This project follows in the tradition of independently implemented transit search algo¬ 
rithms applied to publicly available datasets (such as Petigura et ah 2013b,a; Sanchis-Ojeda 
et ah 2014; Dressing & Charbonneau 2015). These efforts have been hugely successful, 
especially in the field of exoplanet population inference because, thanks to their relative 
simplicity, the efficiency and behavior of these pipelines can be quantified empirically. The 
work described in this Article is built on many of the same principles as the previous projects 
developed for studying Kepler data but our main intellectual contribution is a computation¬ 
ally tractable framework for simultaneously fitting for the trends and the transit signal even 
when searching for planets. 

The Article is organized as follows. In Section 2, we describe our method of extracting 
aperture photometry from the calibrated K2 postage stamp time series. In Section 3 (with 
details in Appendix A), we describe our data-driven model for the systematic trends in the 
photometric light curves and our method for fitting this model simultaneously with a transit 
signal. In Section 4, we give the detailed procedure that we use for discovering and vetting 
planet candidates. To quantify the performance and detection efficiency of our pipeline, we 
test (in Section 5) the recovery of synthetic transit signals, spanning a large range of physical 
parameters, injected into real K2 light curves. Finally, in Section 6, we present a catalog of 
36 planet candidates orbiting 31 stars from the publicly available K2 Campaign 1 dataset. 


'^http://tess.gsf c.nasa.gov/documents/TESS_FactSheet_0ct2014.pdf 
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2. Photometry and Eigen Light Curves 

The starting point for analysis is the raw pixel data. We download the full set of 21,703 
target pixel hies for K2’s Campaign 1 from MAST^. We extract photometry using hxed, 
approximately circular, binary apertures of varying sizes centered on the predicted location 
of the target star based on the world coordinate system. For each target, we use a set of 
apertures ranging in radius from 1 to 5 pixels (in steps of 0.5 pixels). Following Vanderburg 
& Johnson (2014), we choose the aperture size with the minimum CDPP (Christiansen et ah 
2012) with a 6 hour window.^ 

All previous methods for analyzing K2 data involve some sort of “correction” or “de¬ 
trending” step based on measurements of the pointing of the spacecraft (Vanderburg & 
Johnson 2014; Aigrain et ah 2015; Crossheld et ah 2015). In our analysis, we do not do any 
further preprocessing of the light curves because, as we describe in the next Section, we ht 
raw photometric light curves with a model that includes both the trends and the transit 
signal. 

One key realization that is also exploited by the official Kepler pipeline is that the 
systematic trends caused by pointing shifts and other instrumental effects are shared—with 
different signs and weights—by all the stars on the focal plane. For a rigorous theoretical 
analysis of this problem, see (Scholkopf et ah 2015). To capitalize on this, the PDC compo¬ 
nent of the Kepler pipeline removes any trends from the light curves that can be fit using 
a linear combination of a small number of “co-trending basis vectors”. This basis of trends 
was found by running Principal Component Analysis (PCA) on a large set of (filtered) light 
curves and extracting the top few (~ 4) components (Stumpe et ah 2012; Smith et ah 2012). 
Similarly, we ran PCA on the full set of our own generated K2 Campaign 1 light curves to 
determine a basis of representative trends but, unlike PDC, we retain and use a larger num¬ 
ber of these components (150). For clarity, we will refer to our basis as a set of “eigen light 
curves” (ELCs) and the full set is made available online'^. The top ten ELCs for Campaign 1 
are shown in Figure 1. 


^https://archive.stsci.edu/k2/ 

^Note that although we chose a specific aperture for each star, photometry for every aperture radius is 
available online: http: //bbq. dfm. io/ketu. 

^http://bbq.dfm.io/ketu 
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Fig. 1.— The top 10 eigen light curves (ELCs) generated by running principal component 
analysis on all the aperture photometry from Campaign 1. 
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3. Joint transit &: variability model 

The key insight in onr transit search method that sets it apart from most standard 
procedures is that no de-trending is necessary. Instead, we can £t for the noise (or trends) and 
exoplanet signals simultaneously. This is theoretically appealing because it should be more 
sensitive to low signal-to-noise transits and similar methods have been shown to be effective 
for hnding transits in ground-based surveys (Berta et ah 2012). The main motivation for this 
model is that the signal is never precisely orthogonal to the systematics and any de-trending 
will over-£t. This will, in turn, decrease the amplitude of the signal and distort its shape. 
In order to reduce these effects, most de-trending procedures use a very rigid model for the 
systematics. For K2, this rigidity has been implemented by effectively asserting that centroid 
measurements contain all of the information needed to describe the trends (Vanderburg & 
Johnson 2014; Aigrain et ah 2015; Crossheld et ah 2015). In the Kepler pipeline, this is 
implemented by allowing only a small number of PCA components to contribute to the 
fit in the PDC procedure. Instead, we will use a large number of ELCs—a very flexible 
model—and use a simultaneous htting and marginalization to avoid over-htting. 

Physically, the motivation for our model—and the PDC model—is that every star on the 
detector should be affected by the same set of systematic effects. These are caused by things 
like pointing jitter, temperature variations, and other sources of PSF modulation. Each of 
these effects will be imprinted in the light curves of many stars with varying amplitudes 
and signs as a result of the varying flat held and PSF. Therefore, while it is hard to write 
down a physical generative model for the systematics, building a data-driven model might 
be possible. This intuition is also exploited by other methods that model the systematics 
using only empirical centroids (Vanderburg & Johnson 2014; Armstrong et ah 2014; Aigrain 
et ah 2015; Crossheld et ah 2015), but our more hexible model should capture a wider range 
of effects, including those related to PSF and temperature. For example. Figure 2 shows the 
application of our model—with 150 ELCs—to a light curve with no known transit signals 
and the photometric precision is excellent. 

If we were to apply this systematics model alone (without a simultaneous fit of the 
exoplanet transit model) to a light curve with transits, we would be at risk of over-htting 
and decreasing the amplitude of the signal. Figure 3 demonstrates this ehect on a synthetic 
transit injected into the light curve of a typical bright star. The middle two panels in this 
Figure show the light curve de-trended using 10 and 150 ELCs respectively. When only 10 
ELCs are used, the measured transit depth is relatively robust but this model is clearly not 
sufficient for removing the majority of the systematic trends. The model with 150 ELCs 
does an excellent job of removing the systematics but it also distorts the transit shape 
and decreases the measured transit depth, hence reducing the signal strength in the BPS 


spectrum (Kovacs et al. 2002). 


In our pipeline we simultaneously fit for the transit signal and the trends using a rigid 
model for the signal and a relatively flexible model for the systematic noise. Specihcally, we 
model the light curve as being generated by linear combination of 150 ELCs and a “box” 
transit model at a specihc period, phase, and duration. The mathematical details are given 
in Appendix A, bnt in snmmary, since the model is linear, we can analytically compnte the 
likelihood function—conditioned on a specihc period, phase, and dnration—for the depth 
marginalizing out the parameters of the systematics model. The signal-to-noise of this depth 
measurement can then be used as a quality of ht metric or candidate selection scalar. This 
computation is expensive but, as described in the following Sections, it is possible to scale 
the method to a K2-size dataset. The bottom panel of Figure 3 shows the application of 
this joint transit-systematics model to the synthetic transit discussed previously. When the 
joint model is used, the correct transit depth is measnred—the transit is not distorted—bnt 
the systematics are also well-described by the model. 

It is worth noting that this model can be eqnivalently thonght of as a (compntationally 
expensive) generalization of the “Box Least Sqnares” (BLS; Kovacs et al. 2002) method to a 
more sophisticated description of the noise and systematics. Therefore, any existing search 
pipeline based on BLS could, in theory, use this model as a drop-in replacement, although 
some modihcations might be required for computational tractability. 

The choice to use 150 basis functions is largely arbitrary and we make no claims of 
optimality. This valne was chosen as a trade-off between the compntational cost of the 
search—the cost scales as the third power of the size of the basis—and the predictive power 
of the model. In some preliminary experiments, we found that using a larger basis did, as 
expected, lead to a marginally higher sensitivity to small transit signals but the gain wasn’t 
sufficient to justify the added cost. 


4. Search pipeline 

In principle, the search for transit signals simply reqnires evalnation of the model de¬ 
scribed above on a hne three-dimensional grid in period, phase, and dnration, and then 
detection of high signihcance peaks in that space. In practice, this is compntationally in¬ 
tractable for any grids of the reqnired size and resolution. Instead, we can compute the 
values on this grid approximately, but at very high precision, using a two-step procednre 
that is much more efficient. 

Specihcally, we must evaluate the likelihood function for the light curve /„ of star n 
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Fig. 2.— A demonstration of the eigen light curve (ELC) ht to the aperture photometry for 
EPIC 201374602. Top: The black points show the aperture photometry and the green line is 
the maximum likelihood linear combination of ELCs. The estimated 6-hour precision of the 
raw photometry is 264 ppm. Bottom: The points show the residuals of the data away from 
the ELC prediction. The 6-hour precision of this light curve is 31 ppm. Note that although 
we show a “de-trended” light curve to give a qualitative understanding of the model, this is 
not a product of the analysis. In this search for transits, the data are only de-trended for the 
purpose of visualization. 
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Fig. 3.— A comparison between de-trending using different numbers of ELCs and a simulta¬ 
neous fit of the systematics and transit model, (a) The raw photometry for EPIC 201374602 
with a synthetic transit injected at 65 days, (b) The black points show the photometry 
de-trended using a linear combination of 10 ELCs and the green line shows the true transit 
model. In this panel and the next, the maximum likelihood transit depth is computed fol¬ 
lowing BLS (Kovacs et ah 2002). While some of the systematics are removed by this model, 
there is still a lot of residual noise, (c) The same plot as panel (b) but using 150 ELCs to 
de-trend. This model removes the majority of the systematics but also distorts the transit 
and weakens the signal; it reduces the measured transit depth, (d) The final panel shows 
the results of simultaneously fitting for the transit and the systematics using 150 ELCs. The 
maximum likelihood depth (marginalized over the ELC weights) is computed as described 
in Appendix A. Like panel (c), this model removes most of the systematics but does not 
distort the transit or reduce the measured transit depth. 
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given a period P, reference transit time T°, dnration D, and depth Z 

p(f„ I p, r», D, z) . ( 1 ) 

We make the simplifying assnmption that each transit enters this qnantity independently. 
This is not trne; as we change beliefs abont each transit, we change beliefs abont the system- 
atics model, which in tnrn affects the other transits. However, this simplifying assnmption 
is approximately satished for all bnt the shortest periods and leads to a hnge compntational 
advantage. Under this assnmption, this likelihood fnnction can be rewritten as 

p{fjP,T\D,Z) = J] p(/jT^(P, TO), D, Z) (2) 

m=l 

where Tm{P, T^) is the time of the m-th transit given the period P and reference time 
TO, and M{P, T^) is the total number of transits in the dataset for the given P and T^. 
Equation (2) can be efficiently computed for many periods and phases if we hrst compute a 
set of likelihood functions for single transits on a grid in Ti and duration Dk 

{p(fn\T,, D,, ■ ( 3 ) 

Then, we can use these results as a look-up table—with nearest-neighbor interpolation—to 
approximately evaluate the full likelihood in Equation (1). 

In the remainder of this Section, we give more details about each step of the search 
procedure. In summary, it breaks into three main steps: linear search, periodic search, and 
vetting. In the linear search step, we evaluate the likelihood function in Equation (3) on 
a two-dimensional grid, coarse in transit duration Dk and hne in transit time Tm- Then in 
the periodic search step, we use this two-dimensional grid to approximately evaluate the 
likelihood (Equation 2) for a three-dimensional grid of periodic signals. Then, we run a peak 
detection algorithm on this grid that is robust to signals with substantially varying transit 
depths. These transit candidates are then passed along for machine and human vetting. 


Linear search The linear search requires hypothesizing a set of transit signals on a two- 
dimensional grid in transit time and duration. For each point in the grid, we use the model 
described in Section 3 to evaluate the likelihood function for the transit depth at that time and 
duration. Since the model is linear and the uncertainties are assumed Gaussian, the likelihood 
function for the depth (marginalized over the model of the systematics) is a Gaussian with 
analytic amplitude L, mean Z, and variance all derived and given in Appendix A. In 
the linear search, we save these three numbers on a two-dimensional grid of transit times T; 
and durations Dk- The transit time grid spans the full length of Gampaign 1 with half hour 
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spacing and we choose to only test three durations: 1.2, 2.4, and 4.8 hours. Figure 4 shows 
the maximum likelihood transit depth Z as a function of transit time T for the light curve 
of EPIC 201613023, a transiting planet candidate with a period of 8.3 days. 


Periodic search In the period search step, the table of likelihood functions generated in 
the linear search step are used to compute the likelihood of the periodic model (Equation 2) 
on a three dimensional grid in period P, reference time T°, and duration V. At each point 
in this grid, the likelihood function for each transit depth is chosen as the nearest point 
(without interpolation) calculated in the linear search. If the time spacing of the linear search 
is sufficiently hne, this will give a good approximation of the correct periodic likelihood. For 
each periodic model, we compute the likelihood of a model where the transit depth varies 
between transits and the “correct” simpler model where the transit depth is constant. The 
variable depth likelihood is given by the product of amplitudes from the initial search 

M(P,T°) 

p^Mn\P.T^.D) = J] . (4) 

m=l 

Since the likelihood function for the depth at each transit time is known and Gaussian, the 
likelihood function for the depth under the periodic model can also be computed analytically; 
it is a product of Gaussians which itself is a Gaussian 


Pconst(/n \P,T^,D) 


M{P,T°) 

n 

m=l 


2 71 6Zr, 


exp 


V 2SZj J 


where the maximum likelihood depth, for the periodic model, is 


( 5 ) 


M(p,rO) - 

V 2 

Z = az 2^ 




and the uncertainty is given by 


M{P,T°) 


(yz^ 




.1 iZr, 


( 6 ) 


(7) 


Note that this result has been marginalized over the parameters of the systematics model. 
Therefore, this estimate of the uncertainty on the depth takes any uncertainty that we have 
about the systematics into account. 


In general, the variable depth model will always get a higher likelihood because it 
is more flexible. Therefore, a formal model comparison is required to compete these two 
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Fig. 4.— The maximum likelihood transit depth as a function of transit time as computed 
in the linear search of the light curve of EPIC 201613023. After the periodic search and 
vetting this target is fonnd to have a planet candidate with a period of 8.3 days. The first 
transit occurs at 7.4 days on this plot. 
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models against each other on eqnal footing. For compntational simplicity and speed, we nse 
the Bayesian Information Criterion (BIC). The traditional dehnition of the BIC is 

-tsic = lnp(/JF, r», B)-:|liiiV (8) 

where the likelihood fnnction is evalnated at the maximnm, K is an estimate of the model 
complexity and N is the effective sample size. To emphasize that K and N are tnning 
parameters of the method, we rewrite this eqnation as 

1 1 CV 

--BIC = Inp(/JP,T°,Z1)-— (9) 

where J is the nnmber of allowed depths—one for the constant depth model and the nnmber 
of transits for the variable depth model—and a is chosen henristically. For the K2 Cam¬ 
paign 1 dataset, we hnd that a ~ 1240 leads to reliable recovery of injected signals while 
still being fairly insensitive to false signals. 

To limit memory consnmption, in the periodic search, we prohle (or maximize) over T° 
and D snbject to the constraint that BICconst < BICvar and reqniring that the signal have 
at least two observed transits. This yields a one-dimensional spectrnm of the signal-to-noise 
of the depth measnrement as a fnnction of period nsing Eqnations (6) and (7) to compnte 
Z/az at each period. The resnlt is a generalization of the BLS freqnency spectrnm (Kovacs 
et ah 2002) to a light cnrve model that inclndes both a transit and the trends. For example, 
Fignre 5 shows the spectrnm for a planet candidate transiting EPIC 201613023. 

After selecting the best candidate based on the signal-to-noise of the depth, we mask 
ont the sections of the linear search corresponding to these transits and iterate the periodic 
search. This permits ns to hnd second transiting planets in light curves in which we have 
already found a more prominent signal. Under our assumption of independent transits, this 
masking procedure is equivalent to removing the sections of data that have a transit caused 
by the exoplanet that produces the highest peak. For the purposes of this Article, we iterate 
the periodic search until we hnd three peaks for each light curve. This will necessarily miss 
the smallest and longest period planets in systems with more than three transiting planets 
but given the conservative vetting in the next Section, three peaks are sufficient to discover 
all the high signal-to-noise transits. 


Initial candidate list The periodic search procedure returned three signals per target 
so this gave an initial list of 65,109 candidates. The vast majority of these signals are not 
induced by a transiting planet: there are many false positives. Therefore to reduce the search 
space, we estimate the signal-to-noise of each candidate by comparing the peak height to a 
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Fig. 5.— The signal-to-noise spectrum as a function of period for the light curve of EPIC 
201613023. This is the generalization of the BLS spectrum (Kovacs et ah 2002) to this 
simultaneous model of the transit and the systematic trends. To compute this spectrum, the 
results of the linear search (Figure 4) were used as described in Section 4. The top peak (at 
a period of 8.3 days) is indicated with a green dot. Iterating the periodic search found no 
other transit signals above the signal-to-noise threshold. 
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robust estimate of variance in BIC values across period. This is not the same criterion used to 
select the initial three peaks but we find that it produces a more complete and pure sample. 
A cut in this quantity can reject most variable stars and low signal-to-noise candidates that 
can’t be reliably recovered from the data. To minimize contamination from false alarms but 
maximize our sensitivity, we choose a threshold of 15. In absolute value, this threshold is 
somewhat higher than the standard signal-to-noise threshold used when searching for transits 
in the Kepler light curves (for example Petigura et ah 2013a) but given the larger amplitude 
of the systematic noise in the K2 light curves, it is not surprising that a higher threshold is 
required to produce a manageable list of candidates for hand vetting. That being said, it is 
likely that a reduction in this threshold would yield more discoveries at the cost of a larger 
set of hand classifications. 

We also find that the signals with periods < 4 days are strongly contaminated by false 
alarms. This might be because of the fact that our independence assumption (Equation 2) 
breaks down at these short periods. Therefore, we discard all signals with periods shorter 
than 4 days, acknowledging this will cause us to miss some planets (Sanchis-Ojeda et ah 
2014). After these cuts, 741 candidates remain; we examine these signals by hand. The full 
list of peaks and their relevant meta data is available online at^. 


Hand vetting After our initial cuts on the candidate list, the majority of signals are still 
false alarms mostly due to variable stars or single outlying data points. It should be possible 
to construct a more robust machine vetting algorithm that discards these samples without 
missing real transits but for the purposes of this Article, we simply inspect the light curve 
for each of the 741 candidates by hand to discard signals that are not convincing transits. 
The results of this vetting can be seen online®. 

Although de-trended light curves are never used in the automated analysis of the data, 
when conditioned on a specific set of transit parameters, the model produces an estimate of 
what the light curve would look like in the absence of systematic effects. This prediction is 
one of the plots that we examine when vetting candidates by hand. For example. Figure 6 
shows the maximum likelihood light curve for EPIC 201613023 evaluated at the candidate 
period, phase, duration, and depth. Similarly, Figure 7 shows the same prediction folded on 
the 8.3 day period of this candidate. 

After visually inspecting 741 signals, 101 candidate transits pass and are selected as 


®littp: //bbq. dfm. io/ketu 
®http://bbq.dfm.io/ketu 
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astrophysical events. Many of these signals are due to “false positives” such as eclipsing 
binary systems, either as the target star or as a background “blend.” We address this effect 
in the following Section, where we separate the list of candidates into a list of astrophysical 
false positives and planet candidates. 

Astrophysical false positives A major problem with any transit search is the potential 
confusion between transiting planets and stellar eclipsing binaries (EBs). Of particular 
concern are grazing stellar eclipses or stellar eclipses that contribute only a small fraction of 
the total light in a photometric aperture, resulting in greatly diluted eclipse depths able to 
mimic the signals of small planets. 

Ground-based transit surveys have experienced false-positive rates well over 50 percent. 
For example, Latham et ah (2009) reported eight eclipsing binaries and one transiting planet 
among the sample of transit candidates in one held of the Hungarian Automated Telescope 
Network transit search. In fact, the follow-up process to try to rule out such astrophysical 
false positives is a large portion of the effort that goes into a transit survey (e.g., O’Donovan 
et al. 2006; Almenara et ah 2009; Poleski et al. 2010). 

Despite this large fraction of astrophysical false positives in ground-based surveys, the 
primary Kepler Mission saw a much lower false positive rate of only 5-10% (Morton & 
Johnson 2011; Fressin et ah 2013), primarily due to three major factors. First, the superior 
precision of the Kepler photometry enables detection of secondary stellar eclipses, odd-even 
transit depth variations, and ellipsoidal variations (Batalha et al. 2010) to a much lower level 
than ground-based surveys. Second, the relatively small pixels and stable pointing of the 
Kepler telescope has enabled the identihcation of many spatially distinct blended eclipsing 
binaries by means of detailed pixel-level analysis (Bryson et al. 2013) to identify shifts in 
the center of light during transits. And hnally, Kepler is sensitive to much smaller planets 
than ground-based surveys, and small planets are much more common than the Jupiter¬ 
sized planets able to be detected from the ground. We note that while Santerne et al. (2012) 
reported a ~35% observational false positive rate, that study was exclusively focused on 
short-period, large candidates, among which false positives are expected to be more likely. 
Desert et ah (2015) has observationally conhrmed a low false positive rate for the majority 
of Kepler candidate parameter space. 

In K2, the precision of the photometric tests used to vet for such false positives is 
lower and they must be applied with care. There are typically only a handful of transits, 
meaning differences between “odd” and “even” transits must be large to create a signihcant 
difference. Searching for ellipsoidal variations is hindered by the short time baseline and 
the increased photometric uncertainty in K2 data. Centroid variations are feasible in K2 
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time [BJD - 2456808] 


Fig. 6.— The maximum likelihood “de-trended” light curve for EPIC 201613023 evaluated 
at the planet candidate’s period, phase, duration, and depth. The transit times are indicated 
by the green ticks below the light curve. This Figure is only generated for qualitative hand 
vetting and in the search procedure, the model is always marginalized over any choices about 
the systematic trends. 
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but must be treated differently than in the original Kepler mission where this effect was 
generally measured using difference imaging (Batalha et ah 2010; Bryson et ah 2013). 

To do hrst-pass vetting for blended EBs among our catalog of planetary candidates, 
we test for signihcant centroid offsets using the machinery that we have already established 
for modeling the systematic trends in the data, inspired by the methods used to vet Kepler 
candidates (Bryson et ah 2013). If any of the candidates have substantial centroid offsets in 
phase with their transits, this indicates that the signal is likely caused by a background or 
foreground transit of an eclipsing binary and we, therefore, remove it from the hnal candidate 
list. This is only an initial vetting step and a more complete characterization of our catalog’s 
reliability is forthcoming (Montet, etal. in preparation). 

To measure centroid offsets, we start by empirically measuring the pixel centroid time 
series for each candidate by modeling the pixels near the peak as a two-dimensional quadratic 
and Ending the maximum at each time. This method has been shown to produce higher 
precision centroid measurements than center-of-light estimates (Vakili etai, in preparation). 
Figure 8 shows the measured x and y pixel coordinate traces for EPIC 201613023. Much 
like the photometry, this signal is dominated by the rigid body motion of the spacecraft 
and we can, in fact, model it identically. In our analysis, we model the light curve as a 
linear combination of ELCs and a simple box transit model at a given period, phase, and 
duration (Equation Al). Under this model, the maximum likelihood depth can be computed 
analytically. If we apply exactly the same model to the centroid trace, the “depth” that we 
compute becomes the centroid motion in transit in units of pixels. Since the motions won’t 
necessarily point in a consistent direction across transits, we treat each transit independently 
and report the average ofiset amplitude weighted by the precision of each measurement. To 
compute the significance of a centroid offset, we bootstrap the offset amplitude for models 
at the same period and duration but randomly oriented phases. If the centroid measured 
for the candidate transit is substantially larger than the random realizations, we label the 
candidate as a false positive. In practice, the precision of the centroid measurements isn’t 
sufficient to robustly reject many candidates, but two candidates—EPIC 201202105 and 
EPIC 201632708—have offsets 3-a above the median out-of-transit offset amplitude so they 
are removed from the final catalog. For example. Figure 9 shows the in-transit centroid 
offset measured for EPIC 201202105 and compares it to the distribution of out-of-transit 
offset amplitudes. 

A quick a priori estimate of the background blended eclipsing binary rate serves as a 
good sanity check. A query to the TRILEGAL (TRIdimensional modeL of thE GALaxy, 
Girardi et al. 2005) galaxy line-of-sight simulation software reveals that the typical density 
of field stars along the line of sight to the Campaign 1 field is about 7.8 x 10“^ arcsec“^. 
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This gives a probability of about 0.16 that a background star might be blended within a 
8 arcsec radius (two pixels) from a target star. Allowing that ~10% of stars might host close 
binary companions within the period range accessible by this survey, this gives a probability 
of 0.016 that a blended binary star might be chance-aligned within two pixels of any given 
target star. Noting that the average number of planets per star with periods less than 30 
days is about 0.25 (Fressin et al. 2013), we can roughly estimate that we expect <10% of 
our candidates to be caused by nearby contaminating EBs. This estimate suggests that such 
astrophysical false positives should be rare in our sample, consistent with our detection of 
only 2 candidates with clear centroid offsets. 


5. Performance 

To test the performance and detection efficiency of our method, we conducted a suite of 
injection and recovery tests, five per star for all 21,703 target stars. For each test, we inject 
the signal from a realistic planetary system into the raw aperture photometry of a random 
target and run the resulting injected light curve through the full pipeline (except the manual 
vetting). If the search returns a planet candidate—passing all of the same cuts as we apply in 
the main search (except the manual vetting)—with period and reference transit time within 
6 honrs of the injected signal, we connt that injection as necovered. The detection efficiency 
of the search is given approximately by the fraction of recovered injections as a fnnction of 
the relevant parameters. 

To generate the synthetic signals, we use the following procednre: 

1. Draw the nnmber of transiting planets based on the observed mnltiplicity distribntion 
of KOIs (Bnrke et al. 2014). 

2. Sample—from the distribntions listed in Table 1— limb darkening parameters and, for 
each planet, an orbital period, phase, radius ratio, impact parameter, eccentricity, and 
argument of periapsis. 

3. Based on the chosen physical parameters, simulate the light curve, taking limb dark¬ 
ening and integration time into acconnt (Mandel & Agol 2002; Kipping 2010), and 
mnltiply it into the raw aperture photometry. 

We then process these light curves using exactly the pipeline that we use for the light curves 
without injections. Finally, we test for recovery after applying the cuts in signal-to-noise and 
period. We should, of course, also vet the results of the injection tests by hand to ensnre 
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Fig. 7.— The maximum likelihood prediction for the light curve of EPIC 201613023 (see also 
Figure 6) folded on the 8.3 day period of this planet candidate. The points are color-coded 
by time and the median a posteriori transit model is overplotted as a black line. 



10 20 30 40 50 60 70 80 


time [BID - 2456808] 


X 



9.20 


9.15 


9.10 


9.05 


7.1 7.2 7.3 7.4 7.5 7.6 
X [pix] 


Fig. 8.— The centroid motion for EPIC 201613023. Left: The measured x and y pixel 
coordinates as a function of time. Right: The pixel coordinates color-coded by time. As 
identihed by Vanderburg & Johnson (2014), the centroid motions fall in a slowly time variable 
locus. If the centroid coordinates in transit are inconsistent with the out-of-transit motions, 
the candidate is likely to be an astrophysical false positive. 
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Fig. 9.— The estimated in-transit centroid offset for EPIC 201202105 (green line) compared 
to the distribution of 1000 centroid offests computed for randomly assigned phases (black 
histogram). The in-transit measurement is 3-a larger than the median out-of-transit offset 
so it is rejected from the final catalog. 
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that our measurements of detection efficiency aren’t biased by the hand-vetting step but, 
since we chose to limit our sample to very high signal-to-noise candidates, it seems unlikely 
that our hand vetting removed any true transit signals. Any estimates of the false alarm 
rate will, however, be affected by this negligence but we leave a treatment of this for future 
work. 

Figures 10 and 11 show the fraction of recovered signals as a function of the physical 
parameters of the injection, and the magnitude of the star in the Kepler bandpass as re¬ 
ported in the Ecliptic Plane Input Catalog (EPIC^). As expected, the shallower transits at 
longer periods are recovered less robustly and all signals become harder to detect for fainter 
stars. It is worth noting that these Figures are projections (or marginalizations) of a higher 
dimensional measurement of the recovery rate as a function of all of the input parameters. 
For example, this detection efficiency map is conditioned on our assumptions about the ec¬ 
centricity distribution of planets and it is marginalized over the empirical distribution of 
stellar parameters. It is possible to relax this assumption and apply different distributions 
by re-weighting the simulations used to generate this hgure. Therefore, alongside this Ar¬ 
ticle, we publish the full list of injection simulations® to be used for population inference 
(occurrence rate measurements). 

While we argue that the most relevant quantity to use to quantify the performance of 
a transit search pipeline is the efficiency with which it discovers transits, it is also useful to 
consider some other standard metrics. In particular, while de-trended light curves are never 
used at any stage of the analysis, our method does make a prediction for the systematics 
model and we can measure the relative precision of the residuals away from this model. 
These residuals are what would be used as de-trended light curves if that was the goal. 
Figure 12 shows, as a function of the Kepler magnitude reported in the EPIC, the 6-hour 
CDPP (Christiansen et ah 2012) for each light curve after subtracting the best £t linear 
combination of 150 ELCs. 


6. Results 

Out of the 21,703 Campaign 1 light curves, our pipeline returns 741 signals that pass 
the signal-to-noise and period cuts. After hand vetting by the two hrst authors, this list 
is reduced to 101 convincing astrophysical transit candidates. Of these, 36 signals—in 31 


^http://archive.stsci.edu/k2/epic.pdf 
®http://bbq.dfm.io/ketu 
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Fig. 10.— The detection efficiency of the search procedure as a function of the physical 
transit parameters computed empirically by injecting synthetic transit signals into the raw 
light curves and measuring the fraction that are successfully recovered. These tests were 
performed on the entire set of stars so these numbers are marginalized over all the stellar 
properties, including magnitude. 
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Fig. 11.— Like Figure 10, the empirically measured detection efficiency of the search proce¬ 
dure as a function of stellar magnitude as reported by the Ecliptic Plane Input Catalog. This 
never reaches 90 percent because these numbers are marginalized over the range of physical 
parameters shown in Figure 10. Even for the brightest stars, the long period, small transits 
cannot be detected. 
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Fig. 12.— The 6-hour CDPP (Christiansen et al. 2012) for each light curve in Campaign 1 
after subtracting the best fit linear combination of 150 ELCs. For each star, the precision is 
plotted as a function of the Kepler magnitude reported in the Ecliptic Plane Input Catalog. 
The “outliers” in the bottom right corner of the plot are caused by a bright star within the 
photometric aperture and the points in the top left corner of the plot are variable stars where 
the major trends in the light curve are not caused by systematic effects, making the EEC 
model a bad £t. 
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light curves—have no visible secondary eclipse and are deemed planet candidates. These 
planet candidates are listed in Table 2. The two candidates transiting EPIC 201367065 were 
previously published (Crossheld et al. 2015) and the third planet in that system is found as 
the third signal by our pipeline but it falls just below the signal-to-noise cut so it is left out 
of the catalog for consistency. This suggests that a less conservative cut in signal-to-noise 
and more aggressive machine vetting could yield a much more complete catalog at smaller 
radii and longer periods even with the existing dataset. 

The remaining signals are caused by EBs with visible secondary eclipses. In most cases, 
the search reports the secondary eclipse as a candidate and in a few very high signal-to-noise 
cases, the period reported by the pipeline is incorrect and multiple candidates correspond 
to the same transit. It is important to note, however, that the choices made in the search 
were heuristically tuned to hnd planets, not binaries, so our results are not complete or 
exhaustive, especially at short orbital periods. There are other methods specihcally tuned 
to hnd EBs in K2 (such as Armstrong et al. 2014, 2015) and these catalogs contain our full 
sample of EBs and more. 

For the planet candidates, we perform a full physical transit ht to the light curve. To 
do this ht, we use Markov Chain Monte Carlo (MCMC; Foreman-Mackey et al. 2013) to 
sample from the posterior probability for the stellar and planetary parameters taking limb 
darkening and integration time into account. In this ht, we continue to model the trends in 
the data as a linear combination of the 150 ELCs but, at this point, we combine this with a 
realistic light curve model (Mandel & Agol 2002; Kipping 2013a). Even though we have no 
constraints on the stellar parameters, we also sample over a large range in stellar mass and 
radius so that future measurements can be applied by re-weighting the published samples. 
In Table 2 we list the sample quantiles for the observable quantities and the full chains are 
available electronically®. Figure 13 shows the observed distribution of planet candidates in 
the catalog. 

In a follow-up to this Article, we will characterize the stars for each of the candidates 
in detail but for now it’s worth noting that many of the planet candidates are orbiting stars 
selected for K2 as M-type stars. If this rate remains robust after stellar characterization 
and if these numbers are representative of the yields in upcoming K2 Campaigns, the K2 
Mission will substantially increase the number of planets known to transit cool stars. 


®littp: //bbq. dfm. io/ketu 
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Fig. 13.— The a posteriori distribution of planet candidates in the catalog. The error bars 
indicate the 0.16 and 0.84 posterior sample quantiles for the radius ratios. 
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7. Discussion 

We have searched the K2 Campaign 1 data set for exoplanet transit signals. Our search 
is novel because it includes a very flexible systematics model, which is £t simultaneously 
with the exoplanet signals of interest (and marginalized out). By this method, we find 
36 transiting exoplanets, which we have vetted by both automatically and manually and 
characterized by probabilistic modeling. The candidates are listed in Table 2 and posterior 
distributions of planet candidate properties are available^®. 

The flexible systematics model we employ is a 150-parameter linear combination of 
PCA components derived from the full set of 21,703 stellar light curves. That is, it presumes 
that the systematics afflicting each star are shared in some way across other stars while 
the astrophysical signals are causally unconnected (Scholkopf et ah 2015). This assumption 
means that, while there is no formal guarantee that the basis contains no astrophysical 
signals, it is unlikely that the top components will be contaminated. It is our belief—although 
not a strict assumption of our model—that the systematics are caused primarily by pointing 
drifts, or movements of the pixels in the focal plane relative to the stars. In principle, if the 
systematics are dominated by pointing issues, the systematics model could require only three 
parameters—three Euler angles—not 150 amplitudes. However, because (as the pointing 
drifts) each star sees its own unique local patch of flat-field variations, the mapping from 
pointing drifts to brightness variations can be extremely non-linear. Furthermore, because 
when the pointing is moving fast there is a smearing of the point-spread function, there are 
effects keyed to the time derivative of the Euler angles as well. The large number (150) of 
linear coefficients gives the linear model the freedom to model complex non-linear behavior; 
we are trading off parsimony in parameters with the enormous computational advantages of 
maintaining linearity (and therefore also convexity). The computational advantages of the 
linear model are three-fold: Convexity obviates searching in parameter space for alternative 
modes; linear least-squares optimization can be performed with simple linear algebra; given 
Gaussian uncertainties and uninformative priors, marginalizations over the linear parameters 
also reduces to pure linear algebra. 

The goal of this Article was to get exoplanet candidates out of the K2 pixel-level data, 
it was not to generate light curves. That is, both the search phase and the characterization 
phase of the method are approximations to computations of a likelihood function for the pixel 
data telemetered down from the satellite. We did not generate “corrected” or “pre-search 
conditioned” light-curves at any stage; we simultaneously fit systematics and the signals of 
interest to the raw data. For this reason, there is no sense in which this method ever really 


^°http://bbq.dfm.io/ketu 



- 30 - 


produces corrected light curves. 

In this work, we are agnostic about fundamental properties of the host stars. The only 
assumptions we make are that the star targeted by the K2 team is truly the planet host, 
and that there is no dilution by other stars in any aperture. As a result, these posterior 
distributions reflect the maximum possible uncertainty in parameters such as the planet 
radius, which depend sensitively on properties of the host star. To use these distributions 
to characterize the properties of specihc systems, one could re-weight our samples using a 
measurement of the inferred stellar properties. 

This project does not live in isolation and this is certainly not the last time the K2 data 
will be searched! There are other teams searching the K2 light curves for transiting planets 
(A. Vanderburg, private comm.) and they are likely to hnd some planets that we did not 
and vice versa. We make many heuristic choices and short-cuts in this search. For example, 
the choice to work at 150 principal components was based on computational feasibility and 
qualitative tests on a handful of light curves instead of any real model selection or utility 
optimization. 

Another major limitation is that, in principle, the systematics model is designed to de¬ 
scribe spacecraft-induced trends, but not intrinsic stellar variability. In practice, the method 
can still hnd planets around variable stars but a more sophisticated model should be more 
robust in this case. One appealing option would be to model the systematics as a Gaussian 
Process where the input parameters are both time and the same 150 ELCs. Interestingly, 
while this model isn’t linear, the search and marginalization can still be executed efficiently— 
using optimized linear algebra algorithms (Ambikasaran et al. 2014, Foreman-Mackey et al. 
in preparation)—inside the search loop. 

Additionally, while we apply this systematics model simultaneously with a transiting 
planet model to search for planet candidates, this scheme is not restricted to planet searches. 
Any astrophysical event that could be observed in the K2 data could be searched for in the 
same way. By modeling a set of ELCs with any arbitrary data model, events in the K2 
data that appear similar to that data model could be identihed. Such a technique may be 
useful in searching for astrophysical events such as ellipsoidal variations induced by orbiting 
companions, stellar activity, microlensing events, especially in the upcoming Campaign 9, or 
active galactic nuclei variability. 

A substantial caveat to the reliability of all existing transiting exoplanet searches is 
that they all include human intervention. This makes quantifying the false alarm rate of 
these catalogs complicated. There has been some work on automated vetting algorithms 
using supervised classihcation algorithms (McCauliff et al. 2014; Jenkins et al. 2014) but 
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these methods rely on hand classihed examples for training and the performance is not yet 
competitive with hnman classihcation. 

The catalog of planet candidates presented here inclndes only planets with periods longer 
than 4 days and at least two transits in the K2 Campaign 1 footprint. This means that we 
are necessarily missing many planets with orbital periods ontside this range. In particnlar, 
planets with a single transit in the dataset mnst be abnndant. These candidates are the 
most relevant for the stndy of planetary system formation and for statistical inference of 
the distribntion of habitable zone exoplanets. What’s more, given the observing strategy for 
TESS, where each held will only be contiguonsly observed for one month at a time, methods 
for hnding and characterizing planets with a single transit are vital and the new K2 light 
curves are a perfect test bed. 

As a supplement to this Article, we make all the results, data products, and MCMC 
chains available at http://bbq.dfni.io/ketu. The source for this Article, complete 

with the full revision history, is available at http://github.com/dfm/k2-paper and the 
pipeline implementation is available at http: //github. com/df m/ketu under the MIT open- 
source software license. This code and a lot of computation time are all that is needed to 
reproduce the Figures in this Article. 
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Facilities: Kepler 


A. Mathematical model 


We model the raw aperture photometry as a linear combination of 150 ELCs and a 
transit model. Formally, this can be written for the light curve of the fc-th star as 


fk = A Wk + noise 


(Al) 


where 

f k = ( fk,l fk,2 ■ ■ ■ fk,N ) 

is the list of aperture fluxes for star k observed at N times 

/ X T 

t =[ ti h ) 

In Equation (Al), the design matrix is given by 


A = 


( X2,l 

Xl,2 X2,2 


2:150,1 1 "^0(^1) \ 

2:150,2 1 m0(t2) 


\ xi^N X 2 ,n ••• 2:150,AT 1 meitj^) / 


(A2) 


(A3) 


(A4) 


where the Xj^n are the basis ELCs—with the index j running over components and the index 
n running over time—and mg{t) is the transit model 


mg{t) 


— 1 if fin transit 
0 otherwise 


(A5) 


parameterized by a period, phase, and transit duration (these parameters are denoted by 6). 

Assuming that the uncertainties on fj^ are Gaussian and constant, the maximum like¬ 
lihood solution for w is 




(A6) 
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and the marginalized likelihood fnnction for the transit depth is a Ganssian with the mean 
given by the last element of w^* and the variance given by the lower-right element of the 
matrix 

5wk^ ^ ak A) ^ (A7) 

where ak is the nncertainty on fj^. The amplitnde of this Gaussian is given by 

A = “P (“2^? '•'’‘■"‘“’‘■O 

evaluated at the maximum likelihood value Wk*- 
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Parameter Units Distribution 


limb darkening parameters qi and q 2 

orbital period P 

reference transit time T° 

radius ratio Rp/Ri, 

impact parameter b 

eccentricity e 

argument of periapsis to 


— g ~ f/(0, 1) 

days InP ~ U(ln0.5, In 70) 

days TO ~ U(0, P) 

— In Rp/R^ ~ U(ln 0.02, In 0.2) 

— 6~f/(0, 1) 

— e ~ Beta(0.867, 3.03) 

— to ^ U (—TT, vr) 


Table 1: The distribution of physical parameters for the injected signals. The eccentricity 
distribution is based on Kipping (2013b) and the limb darkening parameterization is given 
by Kipping (2013a). 
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Table 2: The catalog of planet candidates and their observable properties. These values and 
their uncertainties are derived from MCMC samplings and the numbers are computed as the 
0.16, 0.5, and 0.84 posterior sample quantiles. The coordinates are retrieved directly from 
the EPIC. 






